Flux networks in metabolic graphs* 
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A metabolic model can be represented as bipartite graph comprising linked reaction and metabo- 
lite nodes. Here it is shown how a network of conserved fluxes can be assigned to the edges of such 
a graph by combining the reaction fluxes with a conserved metabolite property such as molecular 
weight. A similar flux network can be constructed by combining the primal and dual solutions 
to the linear programming problem that typically arises in constraint-based modelling. Such con- 
structions may help with the visualisation of flux distributions in complex metabolic networks. The 
analysis also explains the strong correlation observed between metabolite shadow prices (the dual 
linear programming variables) and conserved metabolite properties. The methods were applied to 
recent metabolic models for Escherichia coli, Saccharomyces cerevisiae, and Methanosarcina barkeri. 
Detailed results are reported for E. coli; similar results were found for the other organisms. 

PACS numbers: 87.16.Yc, 87.18.Vf 



ABBREVIATIONS 

CBM: constraint-based modelling. 

CS: complementary slackness (a property of LP solu- 
tion pairs at optimality). 

GAM: growth-associated maintenance (in relation to 
ATP consumption). 

gDW: gram dry weight (referring to biomass). 
LP: linear programming. 

NGAM: non-growth-associated maintenance (in rela- 
tion to ATP consumption). 



A metabolic network comprises a list of biochemical 
reactions and their associated metabolites As such, a 
convenient representation is in terms of a bipartite graph 
containing reaction nodes and metabolite nodes, with 
edges between nodes indicating that a given metabolite is 
involved in a given reaction [2|. A schematic example is 
shown in Fig. Q] The metabolic network can be modelled 
by chemical rate equations, giving the rate of change of 
the metabolite concentrations in terms of the fluxes, or 
velocities, of the associated reactions. It is widely ac- 
cepted though that the metabolism comes to a steady 
state very quickly, so that the metabolite concentrations 
are unchanging in time. This means that a flux bal- 
ance condition holds, and the set of reaction fluxes (the 
'fluxome') can, essentially, be regarded as the metabolic 
phenotype. Determination of the fluxome is therefore the 
focus of considerable theoretical []J , and experimental ef- 
fort [1, 0| • The global properties of such flux sets have 
been investigated [f|. 

When the network is represented as a bipartite graph, 
the fluxome is traditionally associated with the reaction 
nodes. Here we show how fluxes can be associated with 
the edges of such a bipartite graph, by combining the 
reaction fluxes with any metabolite property that is con- 
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served in the majority of reactions, such as molecular 
weight. Moreover, assuming the flux balance condition, 
such an edge-associated flux network is conserved at all 
the reaction and metabolite nodes apart from a hand- 
ful of sources and sinks. Thus the edge-associated fluxes 
resemble, for example, electric currents in a network of 
resistors @. This observation may help with the visual- 
isation of the flow of material in these complex reaction 
networks. If the flux-balance condition does not hold 
(for example away from steady state), then the edge- 
associated flux network can still be constructed provided 
a set of reaction fluxes is available. In such a case though, 
the edge fluxes are not in general conserved at the reac- 
tion nodes. Finally in the case where the set of reaction 
fluxes arises from the solution to a linear optimisation 
or linear programming (LP) problem, such as commonly 
encountered in constraint-based modelling (CBM), then 
an edge-associated yield flux network can be constructed. 

CBM is now a well-established approach for calculat- 
ing candidate sets of reaction fluxes [l| . It has been ap- 
plied to micro-organisms from all three domains of life 
0> B S EH) and recently extended to encompass hu- 
man metabolism • For the growth of micro-organisms 
a commonly used paradigm has emerged in which the 
metabolic network is augmented with a biomass reaction 
consuming the end-points of metabolism in the appropri- 
ate ratios, and with exchange reactions to represent the 
uptake of substrates and the discharge of metabolic by- 
products. Maximising the flux through the biomass re- 
action amounts to maximising the specific growth rate of 
the micro-organism. This approach has been highly suc- 
cessful at predicting the behaviour of micro-organisms 
[H EH EH, an d nas also been applied to problems in 
metabolic engineering 0, EH ■ 

As already indicated CBM typically leads to an LP 
problem for the set of reaction fluxes. Mathematically, 
every LP problem has an associated dual [H, E3 ■ The 
dual variables are known as shadow prices reflecting an 
economic interpretation of the dual problem. In CBM, 
shadow prices were first investigated by Varma and Pals- 
son who showed that they are, essentially, yield coefli- 
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exchange reactions 



FIG. 1: Schematic metabolic network as a bipartite graph. 
Edges go from reaction nodes (□, blue) to metabolite nodes 
(O, red). Reaction fluxes, v a , metabolite shadow prices, m, 
and the stoichiometry matrix, Sia , combine to make the edge- 
associated yield flux network, J; a . It can be shown (see text) 
that the Ji a are conserved at all nodes except a limited num- 
ber of exchange reactions which serve as sources, and the 
biomass demand reaction which serves as a sink. The biomass 
shadow price is unity, and the flux through the biomass de- 
mand reaction is the specific growth rate /x. 

cients [l8|, Il9| . Later, the dual problem was explicitly 
formulated by Burgard, Maranas, and coworkers, for use 
in multi-level optimisation problems [HI, H3] ■ Recently 
we described a thermodynamic interpretation of the dual 
problem (2lj . The aforementioned yield flux network can 
be generated very naturally by combining a candidate set 
of reaction fluxes with the corresponding set of shadow 
prices. The properties of the yield flux network (such as 
flux conservation) follow from the so-called complemen- 
tary slackness (CS) relations. The parallel construction 
of the various flux networks explains the strong corre- 
lation between shadow prices and conserved metabolite 
quantities such as molecular weight. 



I. METHODS 
A. Conserved edge-associated flux networks 

Mathematically, a metabolic network is conveniently 
described by the stoichiometry matrix Si a , giving the 
number of moles of the z-th metabolite consumed or pro- 



duced by the a-th reaction. We make a distinction be- 
tween balanced 'internal' reactions which typically rep- 
resent biochemical transformations or membrane trans- 
port processes, and imbalanced reactions introduced in 
CBM such as the biomass reaction or the exchange re- 
actions. We suppose there are a = 1 . . . R reactions (of 
all types) and i — 1 . . . M metabolites, with typically 
M < R. The convention we adopt is that Sj a is positive 
for products, negative for reactants. The corresponding 
bipartite graph has R reaction nodes and M metabolite 
nodes (Fig. [I}. An edge connects a reaction node to a 
metabolite node if and only if Si a ^ 0. We additionally 
suppose the bipartite graph is directed and adopt the 
convention that all edges start at reaction nodes and end 
on metabolite nodes. 

In terms of the stoichiometry matrix, the chemical rate 
equations are 

dci 

where Cj is the concentration of the z-th metabolite, and 
v a is the flux through the a-th reaction (reaction veloc- 
ity). The reaction fluxes are typically measured in units 
of mol/gDW.hr where gDW means gram dry weight of 
biomass. In steady state, dci/dt = 0, leading to the flux 
balance condition 

Ea S ia v a = 0. (2) 

Let qt be any property of the z-th metabolite which is 
conserved in the internal reactions, for example molecu- 
lar weight, number of atoms, and so on. Then the corre- 
sponding edge-associated flux network is defined by 

Qia — qiSia^a (3) 

(no implied summation). For example if qi is molecular 
weight we generate what we term the mass flux network, 
if qi is the number of atoms of a given element we gener- 
ate an elemental flux network, and so on. From our sign 
convention for the stoichiometry matrix one has Qi a > 
for edges connected to product metabolite nodes, assum- 
ing qi > and v a > 0. Hence, typically, the flux flows 
from reactants to products. 

The flux network Qi a defined in Eq. [3] is conserved at 
all the metabolite nodes since 

Ec ( 3>c<=«(Ea S i^a) = 0, (4) 

using the flux balance condition Eq. [5J Moreover Qi a is 
conserved at all the internal reaction nodes since 

J2i = qiSia)v a = (5) 

where, by definition, J2i qi^ia — expresses the fact 
that qi is conserved in the a-th reaction. However, Qi a 
may not necessarily be conserved at imbalanced reaction 
nodes since J2i Vi^ia is not necessarily zero. Clearly these 
are important nonetheless, since they represent sources 
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and sinks for the Qi a - As discussed in more detail below, 
in a typical metabolic model these imbalanced reaction 
nodes are restricted to a small set of so-called exchange 
and demand reactions (see Fig. |T|) . 

In some sense these edge-associated flux networks, and 
the yield flux network introduced below, are more funda- 
mental than the set of reaction fluxes. By construction 
they are insensitive to either local rescaling of the reac- 
tion stoichiometry (Si a — > 5; Q x r a and v a — > v a /r a ), or 
local coarse-graining of metabolites (Si a — > Si a x r\ and 
Qi ~~ * 1il r i)- Since these rescalings apply locally like a 
gauge transformation, one might say the edge-associated 
flux networks are gauge invariant. 



B. The linear programming problem in 
constraint-based modelling 

As outlined in the introduction, the typical applica- 
tion of CBM leads to a so-called primal LP problem, 
with a corresponding dual LP problem. The details are 
discussed in the following subsections. 

Primal linear •programming problem : The variables 
in the primal LP problem are the reaction fluxes v a . The 
constraints are the flux balance conditions in Eq. [2J aug- 
mented by additional constraints as follows. Firstly, ther- 
modynamic considerations may lead to some of the in- 
ternal reactions being judged to be irreversible in which 
case v a > 0. In principle the reaction fluxes can also 
be 'capped' (u™ m < v a < v™ ax ) or fixed to a prescribed 
value (v a = vjj x ) although in practice this is rarely done. 

Next, exchange / demand reactions represent the up- 
take or discharge of substrates from the environment. 
By convention they are imbalanced reactions of the form 
M; ^ 0, where Mj is a metabolite, so a positive flux 
represents discharge of the corresponding substrate and 
negative flux represents uptake from the environment. 
For these reactions Si a = — 1. There is no distinction 
between demand and exchange reactions, although we 
tend to restrict the phrase 'demand reaction' to the case 
when the flux is expected to be positive. Exchange / 
demand reactions are classified according to the allowed 
flux range, as follows : 

= (closed), 

> (half-closed; uptake prevented), 

> — v™ ax (open to a limited extent), 
unconstrained (fully-open). 

(6) 

In the third of these, vf 1 ** is the maximum specific uptake 
rate of the given substrate. Most of exchange reactions 
are half-closed, since there is a need to prevent arbitrary 
uptake. For certain essential minerals, dissolved gases, 
nutrients, and vitamins, the exchange reactions are fully 
open so that the corresponding substrates can be freely 
taken up or discharged by the organism. In typical appli- 
cations, one or two exchange reactions are also opened to 
a limited extent, representing growth- limiting substrates 
such as carbon/energy sources. 



The biomass reaction consumes the end-points of 
metabolism such as amino acids, nucleotides, lipids, and 
co-factors. In an extension to the usual paradigm, we 
split this into an irreversible biomass producing reaction 
Ylj 6 2 Mi — > B and an open biomass demand reaction 
B 0, where B is a artificial metabolite representing 
biomass. The stoichiometry coefficients bi typically have 
units mol/gDW. The flux through the biomass demand 
reaction is the specific growth rate, fj,, typically with units 
1/hr. Since B only features in these two reactions, the 
irreversibility of the production step ensures fi > 0. 

Finally, the energetic requirements of the organism are 
taken care of by including growth associated maintenance 
(GAM) and (for high accuracy work) non growth asso- 
ciated maintenance (NGAM) reactions. Again in an ex- 
tension to the common paradigm, we account for these 
by including an irreversible energy-generating reaction 
ATP 4 - + H 2 -> ADP 3 - + H+ + HPO^ + £ where £ is 
a artificial metabolite representing the energy that can be 
gained by hydrolysing (one mole of) ATP. For the GAM, 
£ is included amongst the metabolites consumed in the 
biomass producing reaction with the coefficient b$ repre- 
senting the GAM requirement. For the NGAM, we add 
an energy demand reaction £ with a positive lower 
bound for the flux, V£ > v™ m where v™ m represents the 
NGAM requirement. 

This completes the specification of the primal LP prob- 
lem in typical applications of CBM. The constraints on 
v a specify a so-called feasible solution space. The aim is 
to maximise the specific growth rate, /z, whilst remaining 
within the feasible solution space. The introduction of B 
and £ allows us to move the non-trivial flux bounds and 
the target of the optimisation to demand reactions, with 
a corresponding simplification to the dual problem. Nu- 
merically, solutions to this LP problem can be found by a 
straightforward application of LP techniques, for exam- 
ple the simplex algorithm [H03- A MATLAB toolbox 
for solving LP problems in CBM has been released by 
the Palsson group [26]. For the present study, we used 
a bespoke interface to the GNU LP kit (GLPK), which 
provides an efficient implementation of the simplex algo- 
rithm. 

Dual linear programming problem : Now we turn to the 
dual LP problem. The dual variables are shadow prices 
associated with constraints in the primal LP problem. In 
particular the flux-balance conditions in Eq. [5] generate a 
set of shadow prices Wi for the metabolites. The shadow 
prices in the dual problem are then subject to constraints 
that correspond to the variables (reaction fluxes) in the 
primal problem. As mentioned above, the shadow prices 
can be interpreted as yield coefficients [IH, Ojij]. A full 
derivation of the construction rules presented below is 
given in the Appendix. 

For the internal reactions, the constraints on the iTi 
can be written in terms of the derived quantities [2l| 

Ba = J2i Ki^ia- (J) 

These are defined for all reactions (the resemblence to 
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reaction affinity will be discussed shortly). In the dual 
problem, each reversible or irreversible internal reaction 
generates a constraint on the B a according to : 

R J = Q (reversible, v a unconstrained), 

\ < (irreversible, v a > 0). ^ ' 

This result applies only in the two cases indicated. The 
generalisation to more complicated situations, such as 
fixed fluxes and double-bounded fluxes, is given in the 
Appendix. 

The constraints arising from the internal reactions are 
supplemented by additional constraints arising from the 
exchange and demand reactions. Since these reactions 
involve only one metabolite (M, with Si a = — 1), 
the corresponding constraint simplifies to feature only the 
corresponding metabolite shadow price. The constraints 
associated with the exchange / demand reactions are : 

{= (fully open, v a unconstrained), 

> (half-closed, v a > 0, 

or limited, v a > — w- nax ), 
unconstrained (closed, v a — 0). 

The biomass demand reaction is special and the corre- 
sponding constraint is 

7rg = 1 (biomass). (10) 

This arises because it is the flux through the biomass 
demand reaction that is the optimisation target in the 
primal LP problem. It also make sense since by definition 
the biomass yield coefficient for adding more biomass is 
unity. Given that ttq = 1 is dimensionlcss, the units of TTi 
for the other metabolites are typically gDW/mol (inverse 
to the units of bi). Again this makes sense in terms of 
the TTi being yield coefficients. 

The objective function in the dual LP problem is to 
minimise w — TTivf 1 ax , where the sum is over the lim- 
ited exchange reactions only. One can show from the 
strong duality theorem [1 61 ] that the minimum value of 
w is equal to the maximum value of the specific growth 
rate /i, provided both problems have solutions. It is quite 
common that there is only one limited exchange reac- 
tion, representing single-substrate limitation. In this case 
the dual objective can be taken to minimise the shadow 
price of the corresponding metabolite. It follows that 
v max ^Qgg no j. en ^ er c| ua l problem any more, and the 
metabolite shadow prices are independent of growth rate. 
Also, at optimality one has /i = w = TTivf lax . But /i/ff ax 
is the standard definition of the yield coefficient, confirm- 
ing the interpretation of TTi as the yield coefficient for the 
limiting substrate. 

This construction of the dual LP problem extends and 
simplifies the results presented in [2l|. It also essentially 
recovers the results obtained by Burgard and Maranas 
and coworkers [HI, |2(J. Numerically, the dual problem 
can of course be solved directly, however the shadow 
prices are often obtained 'for free' as a by-product of 



the primal LP solution method. This is the case with 
the simplex algorithm for instance [1 61 ] - 

Complementary slackness relations : At optimal- 
ity a number of so-called complementary slackness (CS) 
relations hold, linking the solutions to the primal and 
dual LP problems. These are also derived in the Ap- 
pendix. For irreversible internal reactions the CS relation 
is B a v a = 0. There is no CS relation for reversible in- 
ternal reactions but, since B a = is imposed, it follows 
that at optimality B a v a = for all internal reactions. 
For the exchange reactions the CS relations are TTiV a = 
for half-closed exchange reactions and TTi(v a + v™ ax ) = 
for the limited exchange reactions. There is no CS re- 
lation for fully open exchange reactions but again, since 
TTi = is imposed (apart from biomass), it follows that 
at optimality TTiV a = holds for all exchange reactions 
except for limited exchange reactions operating at the 
lower flux bound (v a — — uf lax ). There is no CS relation 
for the biomass demand reaction since it is assumed fully 
open. 

C. The yield flux network 

We now show how a pair of complementary solutions 
to the above primal and dual LP problems can be used 
to construct a naturally conserved edge-associated flux 
network, similar to those constructed in section ITA1 We 
start by defining the quantities 

Jia — ^iSi a V a (11) 

which typically have units of 1/hr. The Ji a are conserved 
at metabolite nodes since 

E a Jia = ME a SiaV a ) =0 (12) 

follows from the flux balance condition. At optimality the 
Ji a are also conserved at internal reaction nodes since 

J2i J ia = (J2i "KiSia)v a = B a V a = 0, (13) 

from complementary slackness. 

Exchange reaction nodes are linked by a single edge to 
the corresponding metabolite. For these edges, Sic = — 1 
and Ji a = —TTiVa- Complementary slackness therefore 
implies Ji a = unless the exchange reaction happens 
to be operating at the lower flux bound in which case 
Jia = 7TitiJ nax . For the biomass demand reaction one has 
Jia — —fi since v a = [i and 7rg = 1. 

Thus we conclude that at optimality the yield fluxes 
Jia are conserved at all nodes of the bipartite graph, 
except for exchange reaction nodes where the reaction 
happens to be operating at the lower flux bound, which 
act as sources, and the biomass demand reaction node, 
which acts as a sink. We argue this justifies the notion 
that the Ji a constitute a 'yield flux network' indicating 
how material which contributes to growth is transmit- 
ted through the network. If there is only one limited 
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exchange reaction, the conservation law derived above 
implies 7TjW? lax = fi. But at optimality the strong dual- 
ity theorem shows this is true, as we have discussed in 
section II B[ 



correspondence with conventional non-equilibrium ther- 
modynamics. 



D. Comparison between flux networks 



Genome-scale metabolic reconstructions 



The yield flux network Ji a and the flux networks Qi a 
introduced in section II Al are very similar. They share 
the same small subset of reaction nodes which act as 
sources and sinks, namely the exchange and demand re- 
action nodes (Fig. [IJ) . Thus one might expect the pattern 
of fluxes to be similar. Moreover, since Ji a and Qi a are 
constructed from identical stoichiometry coefficients and 
reaction flux sets (see Eq. [3] and Eq. QT]) , this suggests 
that one would expect a strong correlation between the 
shadow prices 71^ and conserved metabolite properties qi. 
This is confirmed by studies of genome-scale metabolic 
reconstructions, discussed in section [Til We cannot pro- 
vide a proof of an exact relationship between the 7Tj and 
Qi, and indeed the correlation is not expected to be per- 
fectly linear since the corresponding source terms are not 
necessarily in exact proportion. For example the yield 
flux network will typically have a single exchange reac- 
tion node as a source node (i e. the one operating at the 
lower flux bound), whereas the mass flux network has 
sources at all the exchange reaction nodes which carry 
a reaction flux (since all metabolites have some molecu- 
lar weight). By the same argument, it is of course not 
surprising that different conserved molecular quantities 
are only approximately linearly correlated, for example, 
molecular weight is only approximately proportional to 
atom count. 



E. Analogies to chemical thermodynamics 

Recently we described a thermodynamic interpretation 
of the dual problem 2l|. For completeness, we sum- 
marise the analogies and differences between the dual LP 
problem and non-equilibrium thermodynamics as applied 
to these networks by Beard and coworkers [22|, HH, [24[ • 
There is obviously an analogy between Eq. [71 and Eq. [8l 
and conventional chemical thermodynamics [25[ , wherein 



■Ki <-> chemical potential, 
B a «-> reaction affinity. 



(14) 



However the analogy fails to be exact since the CS condi- 
tions require, at optimality, B a v a = 0. This means that 
whenever there is a flux through a reaction (v a ^ 0) the 
corresponding 'affinity' vanishes (B a = 0). This stands 
in sharp contrast to conventional non-equilibrium ther- 
modynamics where a flux through a reaction is usually 
associated with a negative reaction affinity (driving force) 
[22l . I23I [2^ | . It shows that the thermodynamic inter- 
pretation of the dual problem cannot be put into exact 



A genome-scale metabolic reconstruction encompasses 
all the biochemical transformations allowed for by en- 
zymes encoded on the genome of the organism of in- 
terest. As such it represents the entire metabolic ca- 
pability of the organism. A growing number of such 
reconstructions are becoming available. The principal 
genome-scale model used in the present study is iAF1260 
for Escherichia coli [101 ] , which is perhaps the most com- 
plete whole-organism model currently available. We have 
also studied iND750 for Saccharomyces cerevisiae Q , and 
iAF692 for the archeal methanogen Methanosarcina bark- 
eri [9j], in addition to an earlier model UR904 for E. coli 
adjusted slightly to account for later literature @, Ht], Hl| • 

The energetic requirements in these models are repre- 
sented by a GAM component in the biomass producing 
reaction and a separate NGAM reaction. For our calcu- 
lations we retain the GAM requirement, but the NGAM 
requirement was turned off for simplicity (v™ in = 0). We 
have checked that this approximation has little influence 
on our results. Exchange reactions are provided for all 
the extracellular metabolites in these models. By default 
they are half-closed, meaning the flux is constrained to 
be non-negative so discharge only is possible. A subset 
of the exchange reactions are made fully open, on a case- 
by-case basis (details available on request) . In our calcu- 
lations one additional exchange reaction was also opened 
to a limited extent representing the limited availability 
of a substrate under single-substrate growth limitation 
conditions. 

For genome-scale problems it is often the case that even 
at optimality the flux distribution is still not uniquely 
constrained, because there may be alternate pathways in 
the metabolism. Mathematically this shows up in the ex- 
istence of alternate optima in the LP problem [29] . This 
is an interesting phenomenon which somewhat compli- 
cates the analysis. If the primal LP problem has alter- 
nate optima, there will be a similar plurality of solutions 
to the dual problem. However for every optimum solu- 
tion to the primal problem, there exists a complemen- 
tary optimal solution to the dual problem. For instance 
the shadow prices generated 'for free' by the simplex al- 
gorithm are automatically complementary to the primal 
solution. It is important to note that the yield flux net- 
work described above must be constructed using comple- 
mentary solution pairs, since the CS relations apply only 
in this case. The genome-scale models discussed below 
exhibit the phenomenon of alternate optima. However 
we have checked rather carefully that we have used rep- 
resentative solution pairs in presenting our results. 
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G. Statistical analysis 

We undertook statistical analysis of the shadow price 
distributions for selected conditions and organisms al- 
though the results are somewhat inconclusive. We at- 
tempted to fit the observed distributions, using maxi- 
mum likelihood estimators, to a log-Normal, a x distri- 
bution, an inverted x distribution, and a general distri- 
bution ~ exp[— (jrf + LoP)/unr i ] which has tunable expo- 
nential asymptotic behaviour for large and small tti (lu 
and P are parameters). However none of these distribu- 
tions could be said to fit the observed distributions, as 
judged by the Kolmogorov-Smirnov test [3(J. Aside from 
attempting to fit the full distribution, we also determined 
whether the observed distributions were compatible with 
power-law asymptotic behaviour at large and small 7Tj, 
using tail estimators of Hill 3l| and Meerschaert-Scheffer 
[32| . Our results gave exponents systematically greater 
than 3 (see for example Fig. [21 upper plot), which is gen- 
erally taken to be an indication of exponential asymptotic 
behaviour rather than power-law behaviour. When we 
apply these tail estimators to the distribution for | Ji a | , 
we recover an exponent value » —3/2 for large magni- 
tudes, shown as the dashed line in figure [^b). 



II. RESULTS AND DISCUSSION 
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We report results for iAF1260 for E. coli [lOj], under 
various conditions. Similar results were obtained for the 
other organisms and models studied. The data used to 
generate figures [2HU has been compiled into an Excel 
spreadsheet, and is given as supplementary material. 

We first discuss the statistical distribution of the 
shadow prices and yield fluxes. Figure HJa) shows that 
the shadow prices in iAF1260 have a broad distribution 
of around three orders of magnitude. We used statistical 
tests described in section ITGl to analyse the distribution, 
however these were rather inconclusive. We have con- 
cluded though that there is unlikely to be any asymp- 
totic power-law behaviour in the distributions. For the 
cases studied, « 80% or more of the metabolites have a 
positive shadow price, « 15% have a zero shadow price 
meaning that the growth rate is unchanged if the pro- 
vision of these metabolites is altered, and « 5% or less 
have a negative shadow price meaning the growth rate 
actually goes down if that metabolite is injected into the 
system. Figure [^b) shows the distribution of the yield 
fluxes Ji a = TTiSi a v a . This distribution does appear to 
show asymptotic power-law behaviour for large magni- 
tudes. It is notable that the exponent appears to be the 
same as has been found for the reaction flux distribu- 
tion This is interesting since the yield flux network 
is gauge invariant in the sense discussed at the end of 
section II Al whereas the reaction fluxes are not gauge in- 
variant. The fact that we observe power-law behaviour 
in the yield flux network therefore strengthens the earlier 
analysis of [5j . 



FIG. 2: Distribution of (a) shadow prices and (b) yield fluxes 
for E. coli (iAF1260) growing on various substrates under 
aerobic conditions. The dashed line in (b) is the power-law 

| Jia | 

Now we turn to the correlation between shadow prices 
and conserved molecular properties. Figure [3] shows the 
shadow price as a function of metabolite formation free 
energy, molecular weight, and total atom count. To 
obtain these plots, metabolite formation free energies 
(where available) and molecular weights are taken from 
[lOj], and the atom count is computed from the atomic 
formulae in The weakest correlation is with (mi- 

nus) the free energy of formation. This is unsurpris- 
ing since the formation free energy is imperfectly con- 
served in reactions. A stronger correlation is found with 
molecular weight and the strongest correlation is with 
atom count. These quantities are conserved since reac- 
tions in these genome-scale models are charge- and mass- 
balanced. The shadow price is more strongly correlated 
with atom count than molecular weight because there is 
less spread in the magnitude of the values for atom count. 
The shadow prices discussed here are 'molar' yield coeffi- 
cients. One can of course define a 'mass' yield coefficient 
by dividing by the molecular weight. The dashed line in 
figure [3fb) shows that the mass yield coefficients are ap- 
proximately constant with a value of « 0.5gDW/g (we 
have drawn back from undertaking a linear regression 
analysis as we believe this would over-interpret the data 
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FIG. 3: Shadow prices for E. coli (iAF1260) growing on glu- 
cose under aerobic conditions as a function of (a) (minus) 
metabolite formation free energy, (b) molecular weight, and 
(c) total atom count. The dashed line in (b) is 7r; = 0.5 x Mw 
corresponding to a mass yield coefficient of 0.5gDW/g. 



and would not add any new insights). 

The use of shadow prices to measure efficiencies in a 
model of the central metabolism of E. coli was pioneered 
by Varma and Palsson [3 . Our calculations extend 
the scope of this analysis to more recent genome-scale 
metabolic models. Figure [2] and figure 0Ja) shows the 
shadow prices for E. coli grown on four different limit- 
ing carbon/energy sources. Despite the spread in growth 
rates from these sources, by and large there is little dif- 
ference in the shadow price distribution. Invoking the 




(a) atom count 




(b) atom count 



FIG. 4: Shadow prices for E. colt (iAF1260) growing (a) on 
various substrates under aerobic conditions, and (b) on glu- 
cose under aerobic and anaerobic conditions. 



efficiency arguments of Varma and Palsson, this plot sug- 
gests that the metabolic network of E. coli has evolved to 
be equally efficient for growing on a variety of substrates. 
This reflects the 'bow-tie' structure of the metabolic net- 
work [H, [ID, [35[, as substrates are first broken down 
to a dozen or so common precursors, before being re- 
assembled into the components required for growth. 

Figure Hlb) shows a significant overall lowering of the 
shadow prices for anaerobic growth compared to aerobic 
growth. This reflects a reduced efficiency of the network, 
as more effort has to go into satisfying the energetic re- 
quirements of the organism in the absence of oxidative 
phosphorylation. Two further calculations support this 
conclusion (data not shown). Firstly, a similar reduction 
in shadow prices is found for growth under aerobic condi- 
tions with the ATP synthase reaction disabled. Secondly, 
if the NGAM energy demand reaction is thrown fully 
open (£ ^ 0) so that the organism can trivially satisfy 
its energy requirements, the shadow prices are practically 
unchanged on going from aerobic to anaerobic conditions. 



III. CONCLUSION AND OUTLOOK 

To summarise, we have examined the problem of con- 
structing conserved flux networks defined on the edges of 
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bipartite metabolic graphs. We find that such networks 
can be generated by combining a conserved metabolic 
property such as molecular weight, with the reaction 
fluxes. A similarly conserved, edge-associated flux net- 
work can be constructed from a natural combination of 
the primal and dual solutions to the LP problem that typ- 
ically arises in CBM. The correspondence between these 
edge-associated flux networks is responsible for the high 
correlation between shadow prices and conserved molecu- 
lar properties. The construction of these networks opens 
the way for further investigations, both of the global 
properties of the flux distribution (Fig. ^ Q, and the 
scaling theory of transport in complex networks @. 
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APPENDIX 

This appendix presents for completeness a derivation of 
the dual LP problem and the CS relations for the typical 
LP problem that arises in CBM. It is based on the de- 
velopment in §6.5 in [l6| (see also I n this approach 
we assign a Lagrange multiplier to each constraint in the 
primal problem, for example the shadow prices are mul- 
tipliers associated with the flux balance constraints of 
Eq. [U Flux constraints are handled by conversion to 
quadratic equalities. Thus, restricting the analysis for 
the time being to the two most common flux constraints, 
we have 



V a > 



= u 

+ v. 



2 

a s 
max 



(A15) 



The first corresponds to irreversible internal reactions 
and half-closed exchange reactions. The second corre- 
sponds to limited exchange reactions. Adopting the La- 
grange multiplier approach, we replace the original con- 
strained linear optimisation problem by the following 
problem in which we seek the unconstrained maximum 
of 



z = fi + J2 ia nS ia v a + J2' a y a (v a - u 2 a ) 



(A16) 



where the first term is the original objective function /z, 
the second term incorporates the flux balance conditions 
with ni being Lagrange multipliers, and the third and 
fourth terms accommodate the quadratic equalities with 
y a being the multipliers. The prime and double prime 
restrict the sums to the respective reactions with a zero 
or non-zero lower flux bound. We rewrite this as 



+ El VaVa + YlL VotVa 

a V°<K - La V<*< 



(A17) 



where B a — ^2 i TTiSi a is introduced as a definition to 
correspond to the main text. 

From Eq. IA17I one condition that Z is an extremum is 
dZ/dv a = 0, implying 



(unbound), 
B a = { -y a (bound), 
-1 (biomass) 



(A18) 



('unbound' for unbounded reactions; 'bound' for reac- 
tions with a zero or non-zero lower flux bound, 'biomass' 
for the biomass demand reaction which we assume is 
unbounded) . We notice from Eq. IA17I that y a > is 
required for Z to be a maximum, otherwise Z — > oo 
as u a — > ±oo. From this and the second of Eq. IA18I 
we deduce that B a < for reactions with a lower flux 
bound. Taken with the first of Eq. lA18l this gives the B a - 
conditions for the internal reactions quoted in the main 
text. For the exchange reactions, one has S, a = — 1 for 
the metabolite involved in the reaction, and Si a = for 
every other metabolite. Hence B a = — tt,. From this, 
and Eq. IA18I and y a > 0, we recover the conditions on 
the exchange reaction metabolite shadow prices quoted 
in the main text. 

The CS relations follow from dZ/du a = 0. This im- 
plies y a u a = 0, and hence B a v a = 0, or B a (v a — f™ ln ) = 
0, for reactions with a zero, or non-zero, lower flux bound 
respectively. Expanding this to the various cases gives 
the CS relations quoted in the text. 

Still following the development in flr| , we can read 
off from Eq. IA17I that the dual objective function is to 
minimise w = Y^Lv^T^ '■ But Ua = —B a = 7T; for 
exchange reactions, hence w = Yl'i i"i^f lax ) as quoted in 
the main text. 

This whole approach can readily be extended to other 
classes of reactions. For example a reaction with a pre- 
scribed flux leads to the corresponding B a being un- 
restricted and a contribution B a v^ x being added to w 
where is the fixed flux value. As a special case of 
this, a reaction which is disabled (v^ x — 0) has the corre- 
sponding B a made unrestricted in the dual problem. As 
another special case, an exchange reaction with a spec- 
ified flux has a contribution — iriV^ x added to w. It fol- 
lows that -Ki = —dfi/dv^ x - This proves the shadow prices 
are yield coefficients since adding an exchange reaction 
Mj ^ with a fixed flux v a = < corresponds to 
adding the metabolite at a rate |i^ x |. 

Let us finally consider the general case v™ m < v a < 
t>m ax . i 1 ]^ actually specifies two constraints on the flux, 
and thus gives rise to two dual variables. One possible 
interpretation of the resulting dual problem is as follows. 
In addition to B a < one has B a < ^iSi a rather 
than a strict equality. A double contribution B a v™ m + 
(^jiTiSia — -B a )f™ ax is added to w. There are also two 
CS relations, namely B a (v a — u™ m ) = and i^iSia — 

Ba)(C X - V a ) - 0. 
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